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ABSTRACT 

We present results from high-resolution hydrodynamical simulations which explore the 
effects of small scale clustering in star-forming regions. A large ensemble of small-N 
clusters with 5 stellar seeds have been modelled and the resulting properties of stars 
and brown dwarfs statistically derived and compared with observational data. 

Close dynamical interactions between the protostars and competitive accretion 
driven by the cloud collapse are shown to produce a distribution of final masses which 
is bimodal, with most of the mass residing in the binary components. When convolved 
with a suitable core mass function, the final distribution of masses resembles the 
observed IMF, both in the stellar and sub-stellar regime. Binaries and single stars 
are found to constitute two kinematically distinct populations, with about half of the 
singles attaining velocities > 2 km s~^, which might deprive low mass star-forming 
regions of their lightest members in a few crossing times. The eccentricity distribution 
of binaries and multiples is found to follow a distribution similar to that of observed 
long period (uncircularized) binaries. 

The results obtained support a mechanism in which a significant fraction of brown 
dwarfs form under similar circumstances as those of normal stars but are ejected 
from the common envelope of unstable multiple systems before their masses exceed 
the hydrogen burning limit. We predict that many close binary stars should have 
wide brown dwarfs companions. Brown dwarfs and, in general, very low mass stars, 
would be rare as pure binary companions. The binary fraction should be a decreasing 
function of primary mass, with low-mass or sub-stellar primaries being scarce. Where 
such binaries exist, they are either expected to be close enough (semi-major axis ~ 
10 AU) to survive strong interactions with more massive binaries or else born in very 
small molecular cloud cores. 

Key words: accretion - stars: formation ~ stars: mass function - stars: low-mass - 
binaries: general - brown dwarfs 



1 INTRODUCTION 

The advent of infrared imaging has in recent years brought 
about the realization that virtually all stars are formed in 
clusters (e.g. Clarke et al. 2000). On a smaller scale, most 
field stars are known to be binaries (e.g. Duquennoy & 
Mayor 1991) or higher order multiples (e.g. Tokovinin & 
Smekhov 2002), and there is convincing evidence that the 
multiplicity among young stars is even larger (Reipurth & 
Zinnecker 1993; Ghez et al. 1993; Kohler & Leinert 1998). 
Even at earlier stages, Reipurth (2000) reports on a binary 
frequency as high as ~ 80 % for sources driving giant Herbig- 
Haro flows. Additionally, the evidences for coevality of young 
binary systems (e.g. White & Ghez 2001) further suggest 
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that clustering at the scales of molecular cloud cores must 
be a common property of star-forming regions. As shown 
by Sterzik & Durisen (1998), dynamical processes occurring 
within these small stellar groups at birth are likely to afi^ect 
the properties of the resulting stellar objects. Similarly, gas 
dynamics are equally important at this stage: circumstellar 
discs are likely to influence the dynamics of close encounters 
(McDonald & Clarke 1995) and differential accretion from 
a common envelope can imprint a large range of dynamical 
masses on the stars involved (Bonnell et al. 1997). 

Additionally, from a theoretical point of view, the frag- 
mentation of a cloud core into a small aggregate of stars 
can help to attenuate the angular momentum problem in 
star formation (Mestel & Spitzer 1956; Bodenheimer 1995). 
As suggested by Larson (2002), the formation of a binary 
or multiple system can redistribute the excess angular mo- 
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mentum of the infalling envelope into the orbital motion of 
the stellar companions. Other mechanisms such as transport 
processes in discs and outflows might not be as efficient. 

An often used argument against the paradigm of small 
N clustering at birth is based on the sub-millimetre results 
of Motte et al. (1998), who showed that the mass function 
of dense cores in Ophiucus is very similar to the stellar IMF. 
This important result has been widely interpreted as favour- 
ing a one to one mapping from core to star, leaving rather 
limited role for multiple star formation. Since direct evi- 
dence for such substructure is hard to address observation- 
ally (small N clusters disintegrate over a few crossing times) , 
it is interesting to explore the consequences that this hypoth- 
esis can have on the initial mass function and the orbital 
parameters of binary stars, among other observablcs. The 
present paper will focus on small aggregates of N=5 proto- 
stars set in an initially uniform background of molecular gas, 
where the potential of the cluster is initially dominated by 
the gaseous component. Accretion and dynamical interac- 
tions take place naturally in this scenario and provide a di- 
versity of mechanisms by means of which close binary stars, 
multiple systems and brown dwarfs can be formed, and their 
properties statistically compared with observational data. 

The structure of the paper is as follows. In Section 2 
the progress in our understanding of the fragmentation of 
gas clouds and the disintegration of multiple systems is re- 
viewed. Section 3 describes the calculations performed, with 
a caveat on the role of circumstellar discs. Section 4 gives 
a detailed account of the evolution of the clusters. In Sec- 
tion 5 the results per given core mass are presented. Sec- 
tion 6 introduces the convolution with a core mass function. 
The properties of the brown dwarfs formed are discussed in 
Section 7. Finally, our conclusions are given in Section 8. 

2 PREVIOUS WORK 

2.1 The disintegration of small N clusters 

The earliest studies of the disintegration of small clusters 
treated only point mass dynamical interactions (van Al- 
bada 1968a,b). Subsequently, a number of authors have 
investigated the pure N-body problem using a variety of 
numerical codes (e.g. Harrington 1974, 1975; Mikkola and 
Valtonen 1986). In general, the outcome of such simula- 
tions is the break up of the cluster after N crossing 
times, the energy for the break up being released by the 
dynamical formation of a bound subsystem, usually a bi- 
nary star containing the two most massive cluster mem- 
bers. As pointed out by McDonald & Clarke (1993), such 
an outcome has immediate consequences for binary pair- 
ing statistics, implying a predilection for binary member- 
ship amongst massive stars and a dynamical bias against 
brown dwarfs as binary companions. Another obvious im- 
plication of this formation mode is that the kinematics of 
young stars are likely to reflect the hardness of the three 
body encounters that ejected them from the cluster: for a 
cluster containing stars of typical mass M which suffer such 
an encounter at radius rdose, the typical ejection velocity 
is ~ (GM/rdosc)^''^ ~ 1.5r^jJ,(,^kms"-^, where rdoo is nor- 
malised to an encounter at 100 A.U and a mass of IMq. 
Sterzik & Durisen (1995) were the first authors to link ejec- 
tions from compact clusters to the generation of Runaway 



T Tauri stars (RATTs), whilst Armitage & Clarke (1997) 
pointed out that the ejection process would be likely to 
prune the circumstellar discs of such stars and thus render 
them less likely to be Classical T Tauri stars. To date, the 
most accurate and statistically complete study of the pure 
N-body disintegration problem is contained in the series of 
papers by Sterzik & Durisen and collaborators (Sterzik & 
Durisen 1995, 1998; Durisen et al. 2001), which analyse the 
decay channels and resulting kinematics and binary proper- 
tics of a large ensemble of Nbody simulations. This set of 
simulations provides an invaluable benchmark against which 
to compare hydrodynamical simulations. 

It is evident that in the case of young stars with discs, 
encounters at separations less than about 100 AU will not 
in any case be dominated by point mass gravity, but will be 
further modified by dissipative interactions with circumstel- 
lar discs. Clarke & Pringle (1991) and McDonald & Clarke 
(1995) modelled the effect of discs by adding parameterised 
drag terms during close encounters and found that the pos- 
sibility of energy loss via the disc weakened the dynamical 
bias in binary pairing characteristics - in other words, lower 
mass stars ended up in binaries with a higher probability, 
especially as binary secondaries, compared with the results 
of pure Nbody simulations. 

In all the above calculations, the mass of each star is as- 
signed at the outset, so that although they yield predictions 
for binarity and stellar kinematics, they can say nothing 
about the form of the IMF. Such pre-assignment of masses 
corresponds to the case that most of the mass of the par- 
ent core is contained in non-linear density perturbations at 
the onset of gravitational instability. The opposite extreme 
would be a model in which gravitationally unstable fluctu- 
ations contain initially only a small fraction of the core's 
mass, with the remainder being distributed throughout the 
core. In this latter case, the gravitationally unstable seeds 
acquire most of their mass through subsequent accretion, 
which is competitive in the sense that all seeds attempt to 
feed from the same reservoir. In such calculations, the fi- 
nal masses of stars are determined by the same dynamical 
processes that create the binary pairs and lead to stellar 
ejections. Such ideas were first explored in the simulations 
of Bonnell et al. (1997), who showed that due to the in- 
equitable nature of the competitive accretion process, a large 
dynamic range of final stellar masses will result even in the 
case that the seeds all start off with equal masses. More 
recently, Reipurth & Clarke (2001) have specifically linked 
the cluster formation scenario to the production of brown 
dwarfs, arguing that the low mass end of the mass spectrum 
was populated by seeds that had lost out in the competitive 
accretion process and were prematurely ejected. 

2.2 Hydrodynamics: The fragmentation of 
molecular cloud cores 

The favoured mechanism for the formation of most binary 
stellar systems is the fragmentation of a collapsing molecu- 
lar cloud. Fragmentation has been studied numerically for ~ 
20 years (Boss & Bodenheimer 1979; Boss 1986; Bonnell et 
al 1991; Nelson & Papaloizou 1993; BomieU 1994; Burkert 
& Bodenheimer 1993; Bate, Bonnell & Price 1995; Truelove 
et al. 1998). These calculations appear to show that it is 
possible to form binaries with similar properties to those 
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that are observed in fragmentation calculations. However, 
these simulations lack predictive power since the results de- 
pend sensitively on the initial conditions, which are poorly 
constrained. Besides, following the calculation significantly 
beyond the point at which fragmentation occurs is extremely 
computationally expensive. Bate & Bonnell (1997) quanti- 
fied how the properties of a binary system are affected by the 
accretion of a small amount of gas from an infalling gaseous 
envelope. They found that the effects depend primarily on 
the specific angular momentum of the gas and the binary's 
mass ratio. Generally, accretion of gas with low specific an- 
gular momentum decreases the mass ratio and separation 
of the binary, while accretion of gas with high specific an- 
gular momentum increases the separation and drives the 
mass ratio towards unity (see also Artymowicz 1983, Bate 
2000). From these results, they predicted that closer bina- 
ries should have mass ratios that are biased toward equal 
masses compared to wider systems, since the gas which falls 
on to a closer system is likely to have more specific angular 
momentum, relative to the binary's, than for a wider one. 

A more direct approach has recently been conducted 
by Bate et al. (2002a,b), in which a gas cloud which is sub- 
ject to a supersonic turbulent velocity field is modelled. The 
divergence-free velocity field rapidly generates a richly non- 
linear density structure in the gas (Ostriker, Stone & Gam- 
mie 2001). Dense cores are continually created and destroyed 
and, eventually, some of them become Jeans unstable and 
collapse (e.g. see Padoan & Nordlund 2002 for details). How- 
ever, these cores possess internal sub-structure and frag- 
ment further until the opacity limit for fragmentation (a 
few Jupiter masses) is reached. Thereafter, the stellar seeds 
grow in mass by accretion and interact strongly with other 
collapsed fragments and protostellar discs in their vicinity. 
This system readily demonstrates the formation of small N 
ensembles in which the sort of behavior described in the last 
section is observed to occur. Currently, however, the accu- 
racy of the statistical properties of the stars that can be 
determined from such a one-off calculation are limited by 
the small number of objects (« 50) formed. Additionally, 
the computational expense involved in this large-scale sim- 
ulations makes it prohibitive to explore different initial con- 
ditions. Modeling the evolution of many small-N ensembles 
allows us to obtain more accurate statistics and to investi- 
gate the dependence of the results on the initial conditions 
imposed. 



3 THE SIMULATIONS 

The calculations reported in this paper have been performed 
using a three-dimensional smoothed particle hydrodynamics 
(SPH) code based on a version originally developed by Benz 
(Benz 1990; Benz et al. 1990). The smoothing lengths of 
particles vary in time and space, such that the number of 
neighbours for each particle remains approximately constant 
at Nneigh = 50. Wc usc thc standard form of artificial viscos- 
ity (Monaghan &: Gingold 1983) with strength parameters 
Ov = 1 and /3v = 2. 

The stellar seeds arc modelled by sink particles (Bate, 
Bonnell & Price 1995), which interact only via gravity. Sink 
particles are characterized by a constant sink radius Rsink 
such that any gas matter which falls into it and is bound to 



the sink particle is accreted. The gravitational acceleration 
between two sink particles is Newtonian for r > Rsink, but 
is smoothed within this radius using spline softening (Benz 
1990). The maximum acceleration occurs at r ~ Raink/4; 
therefore this is the minimum binary separation Rmin that 
can be resolved. The sink radius has been chosen to be 10~^ 
of the initial radius of the cloud, so that Rmin is ~ 1 AU. 
The gas component is modelled using 10'' SPH particles with 
equal masses. 

3.1 Initial conditions 

One hundred simulations of these small N ensembles has 
been performed. Each simulation consists of an isothermal 
spherical gas cloud, with initial constant density. Five stellar 
seeds (sink particles) were randomly placed inside a sphere 
of radius R, equal to one fourth of the initial cloud radius 
Rci. The choice of N=5 is motivated by two reasons: first, 
we want to study the effect of dynamical interactions in un- 
stable multiple systems (therefore N must be greater than 
2) and second, we would like to allow for the possibility of 
binary-binary encounters to take place (thus, at least N=4). 
The precise choice of N=5, rather than 4 or 6 is not partic- 
ularly relevant, since the results are not sensitive to small 
variations in the number of seeds. Results from clusters with 
a larger number of seeds, say N=10,20, will be provided 
in a future paper. Regarding the initial density distribu- 
tion, the choice of a homogeneous density profile has been 
made following the properties of observed cores: e.g. Alves 
et al. (2001) found that the Bok globule B68 could be fit- 
ted remarkably well by a Bonnor-Ebert distribution, which 
is characterised by a flat inner profile; many other obser- 
vations support this same result (see Andre et al. 2000 for 
a review). In addition, theoretical work by Bate (2000) has 
shown that cloud cores with initially constant density pro- 
files are more likely to produce binary systems with proper- 
ties more similar to those that are observed. A steeper ini- 
tial density profile would enhance dramatically the effects 
of competitive accretion, consecuently producing stellar ob- 
jects with too unequal masses (see Section 4 for details on 
how competitive accretion works.) 

The initial mass of the cloud Mci is a scale- free quantity 
and therefore is set to 1 M© by default, so that the final 
stellar mass gives directly the cloud mass fraction that was 
accreted by that object. Obviously, the final real masses will 
be given by the product of that mass fraction times the 
mass of the core, taken randomly from a certain core mass 
probability distribution. 

The gas component comprises 90 % of the initial mass, 
and is initially in virial equilibrium with its own gravita- 
tional potential. This condition determines the radius of the 
cloud to be equal to: 

Rc. = Igm,^ (1) 

where G is the gravitational constant, jj, is the mean molec- 
ular weight (the gas is assumed to consist of pure molecular 
hydrogen and therefore n is set to 2), Rg is the gas constant, 
T is the temperature, and Mg is thc gas mass. 

Correspondingly, the protostars are also set in virial 
equilibrium with the gas. This condition determines the ini- 
tial velocity dispersion of the seeds. Note that virial equilib- 
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rium means that, in most simulations, if only the gravita- 
tional force between seeds is computed, the majority of the 
protostars are bound to each other, i.e., they form initially 
non-hierarchical high-order multiples. The absolute value of 
the binding energy of the pairs, however, is very low, of 
the order of several x 10~^ times the binding energy of the 
resulting tightest pair. The cluster is bounded by a small 
external pressure to prevent the escape of gas particles at 
the initial stages when the gas thermal energy is comparable 
with its gravitational energy. The gas is initially stationary 
and thus the total angular momentum is equal to zero. The 
initial (3 parameter (ratio of the rotational energy of the 
protostars to the potential energy of the gravitational in- 
teraction between the protostars and the gas cloud) of the 
stellar system has a mean value of 0.03, i.e., the initial an- 
gular momentum of the stellar seeds is not negligible but 
certainly low. Thus, we favour models in which the stars 
lack rotational support against the pull of the gas once this 
collapses towards the centre. 

The numerical values of all the cluster fundamental pa- 
rameters axe given in Table 1. 

An efficiency of 100 % is assumed in all calculations: gas 
is allowed to be accreted until all the cloud mass is in the 
form of stars. Although this is not realistic, the lack of firm 
constraints on star-formation efficiency of molecular cloud 
cores has deterred us from switching the gas accretion off 
at any arbitrary time, thus allowing the investigation of this 
process at all stages of the cluster evolution. 

An unstable multiple is likely to remain in the cluster 
once all the gas has been exhausted. Thereafter, the cal- 
culation is followed using a pure N-body algorithm. This 
has been done for approximately 80 % of the simulations. 
The N-body algorithm used works under the same princi- 
ples as the SPH code: adaptive time-steps evolved using a 
Runge-Kutta-Fehlberg integrator which provides a spatial 
resolution limited by the softening radius e ~ Rsink • No 
rcgularizatioii procedure has been applied. 

Each of these simulations required an average of 
« 450 CPU hours on Sun UltraSPARC 1/5/10 workstations 
for the hydrodynamical calculations. The N-body procedure 
needed an average of ~ 350 CPU hours, the standard de- 
viation of this distribution being much larger than in the 
hydrodynamical case due to the inequitable timescales in- 
volved in the random process of dynamical decay. To test 
the dependency of the results on the numerical resolution, 
calculations were also performed with 5 x 10* particles. The 
results did not show any qualitative difference from the 10* 
particles run. 

These simulations do not include magnetic fields nor rar 
diative transfer. An isothermal equation of state (T = 10 K) 
is used throughout the calculations to resemble the approx- 
imate balance between compressional heating and radiative 
cooling characteristic of the low-density regions of molecu- 
lar clouds. No feedback from the stars has been modelled. 
This may be an adequate choice for low-mass star formation 
where strong stellar winds are not expected. 

3.2 Circumstellar discs 

Protostellar discs are an intrinsic byproduct of the star for- 
mation process. Disc formation occurs whenever any non- 
zero angular momentum is present in the initial gas from 



which the stars form. These discs can be quite large and 
can thus affect the dynamics of the stellar interactions (Mc- 
Donald & Clarke 1995; Hall et al. 1996). The calculations 
reported here do not include any discs in the initial con- 
ditions. The reason for this is twofold: first, to place discs 
at the onset of the simulations would require a large num- 
ber of SPH particles, which would subsequently slow down 
significantly the calculations and therefore make the choice 
of a statistical approach virtually impossible. Secondly, the 
initial conditions in the gas cloud do not provide enough 
vorticity to generate large protostellar discs; initial veloc- 
ity perturbations such as turbulent fields would do this but, 
again, would require a computationally prohibitive number 
of gas particles. Global rotational motion in the cloud is 
known to generate vorticity and will be explored elsewhere. 

The role of protostellar discs is to provide close encoun- 
ters with an additional dissipative term which can affect 
the pairing characteristics of binary stars. Clarke & Pringle 
(1991) and McDonald & Clarke (1995) included discs in 
their model and found that the possibility of energy loss 
via disc weakened the dynamical bias in binary pairs - in 
other words, more lower mass stars ended up in binaries, es- 
pecially as binary secondaries, compared with the results of 
pure N-body simulations. On the other hand, the inclusion 
of gaseous discs introduces a mechanism for the formation 
of stellar objects via fragmentation induced by gravitational 
instability (Bonnell 1994; Whitworth et al. 1995; Burkert, 
Bate & Bodenheimer 1997; Bate et al 2002a, b), process that 
could set a natural scale for the initial mean separation of 
some of the interacting fragments smaller than that assumed 
here. Additionally, circumbinary discs can have an impor- 
tant effect on the evolution of the binary orbital elements: 
the binary may transfer orbital angular momentum to the 
material in the disc via gravitational torques (Artymowicz 
et al. 1991; Bate & Bonnell 1997); this loss of angular mo- 
mentum reduces the separation of the binary. 

The effects of circumstellar disc on the dynamics of 
these small N clusters will be explored in the future. 



4 THE EVOLUTION OF SMALL-N CLUSTERS 

4.1 Onset of instabilities, cloud collapse and the 
formation of a binary system 

As each calculation begins, the protostellar seeds start ac- 
creting gas from the common envelope. The initial gas den- 
sity is constant and there are no pressure gradients to oppose 
the cluster self-gravity, which therefore begins to collapse. 
The stellar motions generate small perturbations in the gas 
density, producing local departures from the expected r~^ 
profile. Such perturbations arise from the accretion of gas 
that is in the vicinity of the stars: under-dense regions are 
created and, in turn, replenished by the infall of gas from 
larger distances. This replenishment is dependent not only 
on the gravitational attraction of the individual stars but 
mostly on the overall gravitational potential of the cluster. 

The presence of gas in the cluster has two main effects 
on the evolution of the system's energy budget (see Bonnell 
et al. 1997). First, there is mass loading, where the g£is is 
accreted on to the individual stars, increasing their masses 
without changing their momenta and hence decreasing their 
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Mel (Mo) 


Rci (AU) Mg (Mo) 


M* (Mo) 


T (K) 


a cr* (km s ^) Rgink (AU) 


Rmin (AU) 


d.,i„i (AU) 


1 


8770 0.9 


0.1 


10 


1 0.2 » 10 


» 1 


10^ 



Table 1. Initial parameters of the simulations. The meaning of the columns is as follows Md ; Total mass of the cloud; Rei: Initial 
radius of the gas cloud; Mgi Initial amount of mass in gas; M»: Initial amount of mass in the stellar seeds; T: Temperature; a: Ratio 
of internal to gravitational energy (a value of 1 means that the cloud contains initially one thermal Jeans mass); ct, ; Initial velocity 
dispersion of the stellar seeds; Rgink^ Accretion radius of the sink particles; Rmin^ Radius of maximum gravitational acceleration; it sets 
the spatial resolution of the calculations; d*_ini: Initial mean distance between stellar seeds 



kinetic energies. Second, the presence of gas acts as a drag 
source on the stars. The motion of the seeds through the 
cluster accelerates and condenses the gas behind the stars, 
which subsequently experiences isothermal shocks that ra- 
diate the additional kinetic energy away. More importantly, 
the cloud collapse proceeds isothermally and therefore the 
gas gravitational energy increases enormously (in absolute 
value), dragging in the stars, until enough accretion on to 
the seeds turns the stellar component into the main focus 
of gravitational attraction. These two processes contribute 
to making the cluster more bound, increasing the chances of 
forming very tight binaries. 

Once the collapse of the cloud drives the stars towards 
the cluster centre, where the density is larger, accretion can 
proceed more efficiently. The accretion rate of a star depends 
predominantly on its position in the cluster. Those stars 
closest to the centre have the benefit of the cluster potential 
in attracting gas. The stars that, at the time of collapse, 
stay at the outside of the cluster never accrete appreciable 
amounts of gas. This differential accretion process will be 
decisive in providing a large range of final stellar masses. As 
soon as one object accretes a significant fraction of the gas 
reservoir, it starts to dominate the potential of the cluster 
and, therefore, turns into the main component of an unstable 
multiple stellar system. Several possibilities are thereafter 
allowed by the chaotic nature of these simulations. 

First, two seeds may start close enough so as to soon 
form a binary, therefore enhancing accretion onto the sys- 
tem and out-weighting the effect of competitive accretion 
between the two components. At this early stages, the seed 
masses arc still close to their initial values and therefore, the 
binary mass ratio is close to unity. Since no initial rotation is 
present in the cloud, the specific angular momentum of the 
infalling gas about the binary rotation axis will be in gen- 
eral lower than the specific angular momentum of the binary. 
This makes the binary mass ratio depart slowly from unity 
and the binary semi-major axis progressively decrease (see 
Section 2.2). Subsequent dynamical interactions can change 
significantly the semi-major ax;is of the binary, produce ejec- 
tion of lower- mass seeds and break up the stellar system into 
a hierarchical multiple, but not change the binary mass ra- 
tio appreciably. The evolution of the separation and mass 
ratio of the two inner components of a binary and a triple 
system are shown in Fig. 1. Note how the semi-major axis of 
the binary changes abruptly when an ejection takes place, 
whilst the mass ratio decreases slowly due to the differential 
gas accretion of the two components. 

Second, it may happen that only one object accretes a 
significant amount of gas before any close binary system is 
formed. In this situation, the other seeds are attracted to- 
wards the most massive member and acquire very eccentric 
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Figure 1. Top left: Time evolution of the separation and, bot- 
tom left, the mass ratio 5 of a binary star system. On the right: 
same as left diagrams but for the inner pair of a triple system. 
The time axis is given in Myr and the separation in AU. The 
dot-dashed lines indicate the ejection times of single stars. As in 
these two cases, the binary mass ratio decreases monotonically 
in most of the simulations. Taking into account that most of the 
binaries end up with a mass ratio in the 0.5-1 range, we conclude 
that the values of q are quite insensitive to differential accretion 
and dynamical interaction processes occurring in the cluster. See 
Section 4.1 for details. 



orbits. The high eccentricity enhances the differential accre- 
tion as only at pericentrc is the density high enough to make 
gas accretion efficient. The result of this process is to drive 
the mciss ratio of the binary to small values, although never 
too far from unity (binaries with mass ratios smaller than 
0.2 are rare). 

Of course, most of the simulations performed display 
an accretion history that lies somewhere in between the two 
extreme cases described above. Yet the common outcome is 
that most of the cluster mass ends up in the binary system, 
as it settles in the gas-rich centre very early on in the cluster 
evolution, and that the mass ratio of this binary depends 
substantially on how soon the pairing occurs relative to the 
accretion history of the seeds. 
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In summary, the consequences of gas accretion are, 
among others, to produce an unstable multiple system fre- 
quently dominated by a binary star whose mass ratio has, 
to a significant degree, been set. From this point of the cal- 
culation onwards (with about 70 % of the gas matter left) 
dynamical interactions play a major role. 

4.2 Dynamical interactions: ejections, binary 
hardening and hierarchical multiples 

Dynamical interactions become important when the cluster 
evolution is close to its first free-fall time. The central binary 
has accreted a significant fraction of the available gas and 
begins to dominate the cluster potential. The other three 
seeds are attracted to its vicinity and very close encounters 
soon take place. The effect of these encounters is twofold: 
first, repeated three-body encounters between the central 
binary and the incoming object contribute to harden the 
binary, therefore decreasing abruptly its semi-major axis, by 
a factor ranging from 2 to 10, and increasing its eccentricity. 
The sudden growth of eccentricity can have a small long- 
term effect on the mass ratio evolution of the central binary, 
as explained in Section 4.1. 

Secondly, the incoming object, with no exception a 
lower-mass member of the cluster, is ejected. Ejection means 
that cither this object stabilizes into a hierarchical orbit (at 
a distance nmch larger than the new semi major axis of 
the binary), or it attains an outward speed larger than the 
cluster present escape velocity, therefore leaving the cluster 
forever. In the latter situation, the ejection halts any further 
accretion, therefore setting the final mass of the oscapcr and 
its kinematical properties. An object which remains bound 
to the cluster in a large orbit may also be deprived of further 
gas accretion if the gas is already too centrally-condensed. 

In general, these simulations show that at least one ob- 
ject and sometimes two arc ejected and escape from the clus- 
ter at an early stage, when about 70 % of the gas matter is 
still present. These escapers leave with a very low fraction of 
the cluster mass and, on average, need only a very soft kick 
to escape the potential well. The typical ejection velocities 
of these objects are ~ 1 km s^^. 

On the other hand, it is common that a non-hierarchical 
triple system survives the cluster break-up once all the gas 
has been exhausted. The decay of this unstable multiple is 
a stochastic process that needs to be followed numerically 
using pure N-body algorithms. Sometimes the triple sys- 
tem disintegrates into a hardened binary and a single. The 
single star, frequently more massive than those ejected at 
previous times, receives also a stronger kick as the distance 
of close approach can be as small as ~ 1 AU (i.e. as low 
as it is allowed by the imnicrical resolution imposed). In 
other runs, the system stabilizes into a hierarchical triple. 
To consider when this is the case, the two-step criterion pro- 
posed by Eggleton & Kiseleva (1995) has been used. First, 
the triple system is partitioned into two pairs; this can be 
done in three different ways (e.g. [l-2]-3, [l-3]-2, [2-3]-l). For 
each such configuration, two instantaneously Keplerian or- 
bits can be computed: the inner orbit of two of the bodies 
in their centre of mass (CM) frame, and the outer orbit of 
the third body around the CM of the inner component in 
the CM frame of the whole system. For the system to be a 
instantaneous hierarchical triple, both orbits must be bound 
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Figure 2. Stellar Mass Fraction Probability Distribution 
(SMFPD): distribution of stellar mass per given core mass, i.e., 
the final mass obtained in the simulations assuming that the ini- 
tial mass of the cores is set to 1 M©. The bimodal shape results 

from the inequitable mass accretion between single stars and outer 
companions (low-mass peak), and binary stars or inner pairs of 
multiple systems (high-mass) peak. 



(negative energy) at least in one of the configurations, and 
the ratio of the outer to inner periods must be greater than 
unity . If there is more than one configuration in which both 
orbits are bound, that with the largest ratio of outer to inner 
periods is selected. Second, if at a given time the triple sys- 
tem is found to fulfill the previous criterion, the ratio of the 
outer pericentre distance to the inner apocentre distance is 
computed. If this ratio is greater than a critical value Yg^'": 



in 



1 + 



3.7 

lout 



2.2 



1.4 



1/3 



1 , 1/3 
1 + 9out 



1/3 1/3 , 1 
9in lout + 1 



(2) 



which depends on the inner qin and outer gout mass ratios, 
then the system is said to be stable for a time 10" x the 
triple outer period, where n is larger than 4 in ~ 95 % of 
the cases. This latter criterion guarantees that a substantial 
fraction of the triples included in our results have survival 
times greater than 100 Myr. 

Finally, approximately 15 % of the simulations produce 
hierarchical quadruple systems, the only escaper being in- 
variably ejected at very early times. Quadruples have been 
identified using a generalized version of the criterion de- 
scribed in the previous paragraph. In very few cases (2 %), 
two mutually unbound binaries are produced, plus one sin- 
gle star. One of these simulations is particularly striking as, 
contrary to the rule, the total mass in the two binaries is 
much smaller than that of the single. 
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Table 2. Fraction of binaries (£5), triples (ft) and quadruples (fq), 
out of the total number of systems (singles + binaries + ...); fm 
refers to the fraction of all multiple systems altogether. 



5 STATISTICAL RESULTS PER CORE MASS 
5.1 Mass function 

Table 2 shows the frequency of binary, triple and quadruple 
systems found in these simulations. The frequency / is de- 
fined as the fraction of binaries/triples/quadruples over the 
total number of systems, where a system can be a single, 
a binary, etc... (e.g. if out of five stars, two are in a binary, 
the binary fraction would be |). 47 % of the simulations 
produce one binary star and three single stars. The binary 
mass ratio is in the range 0.5-1 and the three escapers are 
of much lower mass. In a minority of cases only two or one 
stars (36 % and 15 %, respectively) are ejected as singles, 
with the third (and fourth) being retained at large separa- 
tions. Once again, the binary contains most of the mass, and 
the multiple companion(s) are, like the single stars, of very 
low mass. This pattern of mass acquisition is reflected in the 
mass distribution of stars produced for a given parent core 
mass (see Figure 2). This Stellar Mass Fraction Probability 
Distribution (hereafter SMFPD) is bimodal: one peak cor- 
responding to the binary members and the other to the low 
mass singles and wide companions. 



5.2 Binary mass ratio 

Figure 3 depicts the dependence of binary mass ratio q on 
primary mass. All systems have been decomposed in pairs, 
following the hierarchical distribution. Thus, e.g. a triple 
will contribute to the diagram as two binaries: first, the in- 
ner two bodies (diamonds), and second, the pair formed by 
the centre of mass particle of the inner binary and the ter- 
tiary companion (asterisk). Quadruples are represented by 
squares. Figure 3 illustrates two striking trends that can 
readily be related to the processes of gas accretion and dy- 
namical interactions described in Section 4. Firstly, there is 
a relative paucity of binaries with low mass primaries - es- 
sentially there are no binary primaries with masses less than 
a factor of 0.5 the minimum core mass assumed. This is just 
because the fraction of the core mass that goes into the sin- 
gle stars/outer multiple companions is very small and so the 
binary primary contains a large fraction of the mass of its 
parent core. Secondly, there are very few extreme mass ratio 
systems and, where they exist, the low mass component is 
almost always the wide outermost companion. 

The reason for this second trend is rather different from 
the reason that extreme mass ratio systems are also rare 
in pure N-body simulations. Here, the main effect is not 
so much that low mass companions are exchanged out of 
binaries by more massive ones but that, instead, once two 
seeds form a binary, which very often occurs early on in their 
accretion histories {q therefore being not far from unity), the 
combined gravitational field enhances accretion on to both 
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Figure 3. Mass ratio q versus primary mass, for all the binary 
and multiple stars in the simulations. Multiples have been de- 
composed in bound two-body systems following the hierarchical 
configuration (see text for details). Binaries are represented by 
diamonds, triples by asterisks and quadruples by squares. A large 
fraction of the binaries have mass ratios in the 0.5-1 rajige. Only 
multiples systems can have very low mass companions in outer 
wide orbits. 



components, which subsequently accrete similar amounts of 
gas. 

Likewise, the reason that the binary contains most of 
the core mass is not the result of repeated exchange interac- 
tions, but instead that, once two seeds form a binary, they 
constitute the focus for further accretion, and thereby out- 
strip the other cluster members in their mass acquisition. 

Although these results are roughly compatible with the 
mass ratio distribution derived by Mazeh et al. (1992) for 
close G dwarf binaries, they clearly do not agree with the 
meiss ratio distribution for G dwarf binaries as a whole: the 
results of Duquennoy & Mayor (1991) show a rise with de- 
creasing mass ratio, at least for mass ratios as low as ~ 0.2. 

This deficit can be somewhat alleviated if some of the 
outer pairs in multiples systems are added to the distribu- 
tion. As Figure 4 shows, the rise with increasing mass ratio 
characteristic of close binaries (solid and dot-dashed line) 
contrast with the steep rise in the distribution at very low 
mass ratio when all the pairs are included (e.g. dashed line 
for triples) . This finding is in line with the general feature of 
the simulations: low mass companions are rrmch more likely 
to be the outlier in triple systems than the secondaries in the 
central binary. A prediction of this result is that where wide 
low mass companions have been detected, their primaries 
might, on closer examination, turn out to be binaries (Gizis 
et al. 2001). 

As mentioned in Section 3.2, the presence of circumstel- 
lar discs, neglected in the present simulations, might con- 
tribute to increase the fraction of low mass ratio systems by 
acting as a source of energy dissipation in stellar encoun- 
ters, or by fragmenting into one companion when orbiting 
a single star. Future experiments (Delgado-Donate et al. in 
prep.) will focus on larger N ensembles in order to discover 
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Figure 4. Mass ratio q distribution. The thin solid and dot- 
dashed line stand for pure binaries and inner pairs in multiples, re- 
spectively; the dashed line represent outer pairs in triple systems, 
and the double-dot-dashed line outer companions in quadruples. 
The binary fraction is an increasing function of q, except for 
higher-order systems, characterized for extreme values of q. Note 
how the distribution for pure binaries is very shallow (fiat to 
first order) , whilst that of inner pairs of multiples increases more 
steeply, due to the fact that tertiary companions are in average 
less massive than single objects ejected from the cluster - the 
reason for this is that a multiple system can achieve a stable hier- 
archical distribution more easily if the outermost companion has 
a much smaller mass than each of the binary components. 



whether ejections due to binary-binary interactions can give 
additional insight into the formation of extreme mass ratio 
systems. Likewise, the evolution of N=2 systems (studied 
by Bate 2000 for circular systems), in which competitive ac- 
cretion is the dominant process and there is no place for 
dynamical interactions, will be explored. 

5.3 Semi-major axis distribution 

Figure 5 depicts the semi-major axis dependence on primary 
mass. As it can readily be seen, most binaries (diamonds) 
are rather tight, with a mean semi-major axis (a) w 5 AU. 
On the other hand, wide systems (asterisk for triples, and 
squares for quadruples) correspond, with few exceptions, to 
triple and quadruple stars. This result can be immediately 
related to the pattern of dynamical decay outlined in Sec- 
tion 4.2. Given that a multiple system will only survive if its 
orbital elements adopt particular ranges of parameter space 
(so that it is sufficiently hierarchical to be dynamically sta- 
ble; see e.g. Eggleton & Kiseleva 1995), triples and quadru- 
ples must display a semi-major axis distribution whose mean 
is > 10 times larger than that of the binaries comprising the 
inner components. 

The minimum semi-major axis for binaries is set by the 
numerical resolution achieved in these calculations, as pa- 
rameterized by Rmin- On the other hand, the low mean value 
of the semi-major axis distribution is determined by the 



Figure 5. Semi-major axis a versus primary mass. Symbol code 
as in Figure 3. Mostly close binaries are formed in these simula- 
tions: the mean a is smaller than 10 AU, but follows a Gaussian 
probability distribution in logarithmic space (see Duquennoy & 
Mayor 1991). This small (a) is due to the initial centrally concen- 
trated configuration of the stellar seeds and the number of close 
dynamical interactions occurring in the cluster (in turn correlated 
with the number of protostars N=5 used in these models). On the 
contrary, multiples have large a, since the components must follow 
a hierarchical distribution, with flouter /dinner being much greater 
than 10. 



number and strength of those dynamical interactions that 
lead to an ejection. Therefore, (a) is strongly anti-correlated 
with the number of seeds N—5 used throughout these calcu- 
lations. The initial centrally concentrated spatial distribu- 
tion of the protostars (the initial separation between seeds 
is ~ 10^ AU, 10 X smaller than the core radius), is also re- 
sponsible of the formation of very close binaries. A more dis- 
tributed initial configuration would hace produced a larger 
fraction of wide binaries. 

From these results, together with those described in Sec- 
tion 5.2, it is clear that there exists a relation between the 
mass ratio, semi-major axis and multiplicity of the system. 
In this kind of simulations, wide binaries (semi-major axis 
> 10"^ AU) could be formed if there is no orbital decay due to 
dynamical interactions. For this to happen the binary com- 
ponents must accrete mass in relative isolation and there- 
fore, companions at larger separations are not likely to be 
present. Since there is no other mechanism that can reduce 
the orbit substantially (the collapse of the core induces a 
semi-major axis shrinkage of ~ an order of magnitude), the 
two stars remain too distant to compete equally for the gas 
reservoir, consequently driving the mass ratio quickly to very 
low values. 

5.4 Eccentricity 

Among binary stars (i.e. excluding all outer pairs in mul- 
tiples), the eccentricity cam take a whole range of values 
from 0.1 to 1, as shown by Figure 6. The thick solid line 
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Figure 6. Distribution of orbital eccentricities e. The thin solid 
and dashed line correspond to double stars and inner pairs of 
multiple systems, respectively. The dot-dashed line refers to triple 
and quadruple stars. The thick solid line histogram includes all 
the pairs. A linear fit is shown to be a reasonable first-order rep- 
resentation of the total e distribution. 



histogram includes all the pairs, whereas the thin solid line 
refers to pure binaries, the dashed line to the inner binaries 
in multiples, and the dot-dashed line correspond to triple 
and quadruple systems. It is apparent that there is a pref- 
erence for high eccentricity values, following a probability 
distribution with a functional form f(e) oc e. Although the 
eccentricity is the orbital parameter most sensitive to star- 
disc interactions (high eccentricities would be damped by 
the accretion of material from the disc), this result is in 
good agreement with the expected distribution of eccentric- 
ities for close young binaries (iu which tidal circularization 
has not had time to operate yet; Ambaxtsumian 1937), and 
long period binaries (Duquennoy & Mayor 1991). 

Multiple systems arc characterized by larger eccentrici- 
ties, as it is expected for objects which have undergone sev- 
eral close triple encounters before settling into a stable hier- 
archical configuration. This trend is even more accentuated 
in quadruples than in triples. Thus, one would expect a mul- 
tiple system to display an increasingly larger eccentricity as 
one moves along the hierarchy, as well as a progressively 
smaller mass ratio. This pattern is reflected in the diagram 
shown in Figure 6: the eccentricity distribution of inner bi- 
naries (dashed line) is approximately flat whilst that of outer 
pairs (dot-dashed line) peaks at very high eccentricities. 

5.5 Kinematics 

The break-up of these small aggregates imprints a clear kine- 
matic signature on the objects involved. Binary stars attain 
centre of mass velocities that, in average, arc a factor of 
10 smaller than the typical ejection speed of single stars 
(see Figure 7; singles are denoted by crosses and multiples 
as in previous figures; the dashed and dotted lines indicate 



Figure 7. Velocity in km in terms of mass (in Mq). Crosses 
indicate single stars. Binaries and multiples are denoted by the 
same symbol code as in Fig 3; for these, the depicted mass is the 
primary mass of the pair. The dashed and dotted line indicate 

the mean velocity of multiples and singles, respectively. Above 
the dot-daslicd line arc located those objects with a speed larger 
than 4 km s^'^. The velocity offset between multiple systems and 
singles is approximatelly one order of magnitude. There is not a 
clear dependence of ejection velocity on mass for single objects. 



the multiples and singles mean velocity, respectively). This 
distinction between the velocity distribution of singles and 
binaries is in essence due to the fact that single stars are all 
ejected from the core, whereas the binary (or multiple) star 
remains close to the core centre of mass. A notable differ- 
ence between the N-body case (e.g. Sterzik & Durisen 1998) 
and these hydrodynamical simulations is the relationship be- 
tween the typical ejection speeds and the initial parameters 
of the core. In the dissipationless calculations, the final ve- 
locities of ejected staxs are of the same order as the virial 
velocity of the initial cluster. In the simulations reported 
here there are also a few objects which are ejected early 
on with such velocities, where the initial virial velocity of 
the core is ~ 0.2 km s^^. However, most of the stars arc 
ejected with velocities much greater than this: about half 
attain speeds greater than 5 times the initial virial veloc- 
ity (the corresponding number for Sterzik & Durisen is less 
than 1 %), and there is also a significant minority that at- 
tain velocities > 4 km s~^ (dot-dashed line in Figure 7). 
This may be readily understood inasmuch as in gas dynam- 
ical simulations the protostellar system shrinks as the gas 
is accreted, so that later interactions are much closer and 
produce correspondingly larger velocity kicks. 

Another noteworthy feature of the kinematics of ejected 
objects is shown in Figure 8. The velocity of singles has 
been averaged in ten logarithmic bins (squares) spanning 
the whole range of masses (per given core mass). The best 
fit to these statistical averages does not display any signif- 
icant correlation between ejection speed and escaper mass, 
within the error bars (given by Poisson noise). This result 
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Figure 8. Speed in 1cm versus mass, for single stars only. The 
velocities have been averaged in ten logarithmic bins (squares) 
and the best linear fit computed. Within the error bars (stan- 
dard deviations of the distribution inside every bin), there is no 
significant variation of ejection velocity with mass. 



Figure 9. Ejection speed in km versus the mass in the stellar 
component at the time of ejection. The stellar mass is shown in 
units of the initial core mass. A clear correlation between these 
two quantities can be seen, as it is expected if the distribution of 
escaping velocities depends on the binding energy of the cluster 
at the time of ejection. 



arises naturally from the dynamics of the accreting cluster. 
Low mass objects can be ejected at all stages of the cluster 
evolution since they can remain bound to the cluster at large 
distances from the gas-rich centre; thus, their speed pattern 
will reflect the different degrees of impulse needed to over- 
come the cluster potential as it gets stronger. This modus 
operandi is illustrated by the laxge spread in velocities shown 
in Figure 8. Note that the apparently lower dispersion in 
ejection velocities for higher mass singles (mass > 0.06 M©) 
results from the smaller rmmber of objects located in this 
mass interval, as demonstrated by a Kolmogorov-Smirnov 
test applied to the speed distribution. 

Consequently, the escaper's speed depends mainly on 
the binding energy of the cluster at the time of ejection 
(see Bate, Bonnell & Bromm 2003 for a similar conclusion), 
which in turn would correlate roughly with the total mass 
in the stellar component at that time. This correlation can 
be seen in Figure 9. 

It is unclear however whether this predicted difference 
in the velocity of young single and binary stars would in 
practice be detectable in star forming regions, since the ve- 
locity offset between the two (a few km s ~^) is comparable 
with the corc-corc velocity dispersion in such regions. Ad- 
ditionally, the presence of various binaries in the same core 
might spread the velocity dispersion of the binaries via ejec- 
tion of some of the pairs, therefore making it less distinct 
from the population of escapers (Bate, Bonnell & Bromm 
2003). 



6 THE CONVOLUTION 

In order to relate the findings described in Section 5 to the 
resultant IMF, it is necessary to convolve the final stel- 



lar masses per given core mass (or SMFPD) with a core 
mass function (henceforth CMF). This means that the final 
masses of the protostars and related quantities are deter- 
mined by two processes: first, the seeds gain mass according 
to the acquisition pattern found to operate in these cores; 
and second, the star obtains its final mass once the initial 
mass of its parent core is determined. In this operation we 
are assuming that the division of mass is a scale- free process. 
This would certainly be the case either if there were no sys- 
tematic dependence on core mass of dimensionless variables 
such as e.g. the Jeans number (as controlled by equation 
(1)), or else if the results depicted here are insensitive to 
those variables. 

Mathematically, this convolution differs from the usual 
one in the sense that two numbers (the stellar mass fraction 
of a certain seed and the mass of the core in which this star 
is embedded) axe combined by a multiplicative instead of 
additive operation. That is, the final mass of a stellar object 
m, will be given by the product r, x Mc, where r, is the 
mass accreted by the seed in units of a given core mass, and 
Mc is the mass of that core. 

Unfortunately, the choice of a core mass function is not 
unique since the process that leads to the formation of dense 
star-forming cores is only partially understood. Several func- 
tional forms for this CMF are possible: 

(i) A power-law with index of —2, as expected from the 
equipartition of bulk kinetic energy and gravitational energy 
in a pressureless fluid. 

(ii) A log-normal mass function, as derived from the 
density probability distribution function of Jeans unstable 
clumps in computer simulations of compressible magnetohy- 
drodynamic turbulence (Padoan et al 1997; Klessen 2001). 

(iii) A power-law with Salpeter or nearly-Salpeter slope 
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Figure 10. Initial Mass Function (IMF) for the stars and brown 
dwarfs formed in the simulations. This IMF is the result of the 
convolution of the bimodal SMFPD with a core mass function 
(CMF) taken to be a power-law of slope a=— 2.35 (Salpeter) and 
cutoff at 0.25 and 10 Mq. For masses approximately above the 
product of the largest SMFPD value with the lowest CMF value, 
the mass function resembles the observed IMF. Below, the distri- 
bution is approximately flat in logarithmic scale (slope — 1). 



(a = —2.35), as found by Mottc et al. (1998) in the mass 
distribution of dense cores in Ophiucus. 

(iv) A combination of power-law and log-normal distri- 
butions, as proposed by Padoan & Nordlund (2002) based 
on turbulent fragmentation calculations. 

In order to fully understand the effects that a multi- 
plicative convolution has on the CMF, wo must note that a 
SMFPD described by a delta function would conserve the 
shape and slope of the CMF, and simply displace it to the 
position of the delta. If instead of a delta function we use a 
narrow box function, the CMF is additionally stretched to 
the new limits given by [ rt,min x Mc,min , r«,max x Mc.max ]. 
Also, the convolution of two power-laws - e.g. a SMFPD of 
slope —a confined between and 1, and a CMF of slope — /3 
with an arbitrary larger upper cutoff (whore a is assumed 
to be smaller than B) would, in general, give rise to a 
new function characterized by a break at a position given 
by the product of the lowest CMF value and the highest 
SMFPD value. This break would split the new function into 
two regions: at lower values, a power-law of slope —a, and 
at higher values, a power-law of slope — /3. Note that a and 
P would exchange their range of influence if we choose a to 
be larger than (3; that is, the steepest slope dominates at 
larger values. 

6.1 Initial Mass Function 

Figure 10 shows the result of a convolution of the bimodal 
SMFPD with a CMF taken to be a power-law of index -2.35 
(Salpeter IMF) in the range of core masses 0.25 to 10 M0. 
The choice of the power index and upper cutoffs are moti- 



vated by the millimetre continuum observations of Motte et 
al. (1998) and the H"C0+ survey of Onishi et al. (2002), re- 
spectively. The lower cutoff comes from the minimum mass 
a core can have so that the initial seed mass lies above the 
opacity limit for fragmentation (~ 5 Jupiter masses; Low & 
Lynden-Befl 1976; Rees 1976; Boss 1989; Bate 1998, 2002), 
given the initial conditions imposed - 10% of the mass is in 
stellar form at the start of the calculations. Clearly, the re- 
sulting IMF essentially follows the core mass function down 
to a stellar mass close to the minimum core mass. The reason 
for this correspondence is that for a steep core mass func- 
tion, the majority of stars of a given mass belong to the high 
mass peak of low mass cores, rather than the low mass peak 
of the less numerous liighcr mass cores. Thus the rmmber 
of stars of given mass is governed by the frequency of cores 
whose mass exceeds the stellar mass by a rather modest fac- 
tor (less than a factor two) and the IMF therefore follows 
the core mass function. 

Similarly, convolutions performed using log-normal 
mass functions as CMF, in which the characteristic mass 
and standard deviation adopt typical values, and the cutoffs 
are set at the same values chosen for the power-law CMF 
discussed above, produce invariably the same outcome. That 
is, the resulting IMF resembles the original log-normal mass 
function and simply differs in its width, which stretches to 
the values set by the products of minimum and maximum 
core and stellar masses. Also, as it may be expected from the 
previous results, the functional form proposed by Padoan 
& Nordlund (2002) for the CMF (in the original, IMF) is 
equally preserved. 

Thus we find that although star formation in this case is 
far from being a one to one mapping between core and star 
mass, the IMF of the stars is controlled by, and nearly iden- 
tical to, the core mass function. We note that the observed 
similarity between the IMF and core mass function cannot 
therefore bo used as an argument against multiple fragmen- 
tation and competitive accretion inside star-forming cores. 

6.2 Binary Properties 

The convolution of core masses with a CMF has an effect on 
the resulting properties of binary stars. Although the binary 
mass ratio and eccentricity is independent of the individual 
masses of the binary components, the dependence of these 
parameters on the mass of the primary is modified by the 
CMF assumed. Additionally, if the initial mass of a core 
and its radius are assumed to be correlated in some way (as 
it is the case when e.g. the initial cloud conditions are set 
according to equation (1)), the semi-major axis distribution 
is also affected by the choice of the CMF. Equation (1) also 
states that the final velocity dispersion is independent of the 
initial mass assumed for each individual core. 

Assuming (1) to hold, and making use of the same CMF 
as in Section 6.1, wo have derived the scaled primary masses 
and semi-major axis of all the pairs. Fig 11 shows the semi- 
major axis dependence on primary mass. The linear rela- 
tion between primary mass and semi-major axis imposed by 
equation (1) is reflected in the diagram, as well as the effect 
of the convolution, which has smeared out the distribution 
of points initially clustered around the original mean value. 
From this diagram, it can also be seen that the process of 
convolution with a CMF results in a large fraction of bina- 
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Figure 11. Same as in Figure 5, after performing a convolution 
with a power-law CMF. 
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in the model. Systems where repeated strong interactions 
are not so likely might produce longer period binaries - for 
this to be a solution to the lack of wide binaries in these 
simulations, the fraction of cores with two initial seeds must 
be comparable or larger than that of cores possesing more 
substructure. The addition to the gas of some initial angu- 
lar momentum might as well enhance the fraction of binaries 
with longer periods. 

In summary, the choice of an observationally-based 
CMF results in a sot of binary masses and orbital param- 
eters that are in reasonable agreement with observations. 
The formation of wide binaries, as suggested above, might 
be accomplished in a scenario where dynamical interactions 
were not important, as in cores where only two or three col- 
lapsed fragments are present. The origin of very low mass 
ratio binaries, however, remains obscure: both pure N-body 
calculations (e.g. Sterzik & Durisen 1998) and the present 
hydrodynamic- N-body models fail to produce as large a frac- 
tion of extreme mass ratio binaries as observational data 
suggest (the Duquennoy & Mayor mass ratio distribution 
peaks at a value of 0.2). Although some of these double stars 
might turn out to have internal substructure, in the sense 
of the inner body being in fact a binary star (e.g. Tokovinin 
& Smekhov 2002), - this would be immediately accounted 
for in our simulations -, it is obvious that very low mass ra- 
tio binaries constitute a significant population and its origin 
needs to be further investigated. 



\ 10 7 




10^ 10^ 10^ 10"^ 10^ 

PERIOD [Days] 

Figure 12. Distribution of orbital periods for the binary stars 
produced in these simulations. Shown is a Gaussian fit to the 
histogram. 



ries with low mass primaries - binary primaries have a mass 
~ 0.5 X the core mass assumed, but there are many more 
cores with masses below 1 M0 than above. 

The histogram of periods is shown in Fig 12. It can be 
fitted by an approximate Gaussian distribution which peaks 
at ~ 10^ days. Patience et al. (2000) find a similar value 
for binaries in dense environments (open clusters), whereas 
Duquennoy & Mayor (1991) show an observational maxi- 
mum at ~ 10** days for field G type stars (comparable to 
that of T Tauri stars). This discrepancy might be due to 
the number of stellar objects imposed as initial conditions 



7 BROWN DWARFS 

Once a convolution with a suitable CMF is performed (e.g. 
a power-law with nearly Salpeter slope within the observed 
mass range), the calculations show that the population of 
escapers and outer companions is abundant in brown dwarfs. 

As described in Section 4.2, low mass objects are nat- 
urally produced when dynamical interactions within an un- 
stable accreting multiple are taken into account. These in- 
teractions provide a simple dynamical explanation for the 
formation of sub-stellar objects, which start with a mass 
of the order of the opacity limit of fragmentation and are 
prevented from reaching the hydrogen burning limit due to 
the limited amount of gas they can accrete before they are 
ejected from the common envelope (Reipurth & Clarke 2001; 
Bate et al 2002a). 

Figure 13 shows the mass ratio of multiple systems in 
terms of the mass of the primary component. As explained 
in Section 5.2, all multiples have been partitioned in two- 
body systems (symbol code as in previous figures). Below 
the solid line are located those pairs which posses a brown 
dwarf companion; binary brown dwarfs occupy the region to 
the left of the dashed line. It can be readily seen that brown 
dwarfs produced by dynamical interactions are infrequent as 
binary companions. On the other hand, among hierarchical 
triples and quadruples the outermost object is very often a 
brown dwarf, with a probability higher than 90 % for low- 
mass primaries - this fraction is certainly an upper limit, 
since some of the triple and quadruple systems included in 
this study might not be stable enough to survive the pre- 
main-sequence stage. In this sense, no brown dwarf desert at 
wide separations should be apparent if the inner components 
were not resolved (see Gizis et al. 2001). Brown dwarfs can 
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Figure 13. Same as in Figure 3, after applying a convolution 
with a CMF as described in text. Note that here both axis are in 
logajrithmic scale. 



Figure 14. Same as in Figure 7. The masses have been convolved 
using a CMF as described in the text. 



be wide binary companions, but in this case the primary is 
predicted either to bo a binary system or to have a low mass 
(less than fs 0.3 Mq). Binary brown dwarfs are equally im- 
probable in this scenario, with the exception of some scarce, 
tight, nearly equal mass pairs. 

The above results on brown dwarf binary pairings rep- 
resent a critical area for confronting dynamical decay mod- 
els with observations. It is currently unclear how the binary 
fraction changes in the case of very low mass primaries, since 
the detection of brown dwarf binaries is in its infancy (sec, 
for example, Reid et al. 2001). If the binary fraction were 
indeed lower for brown dwarfs, it could readily bo explained 
by such models, whereas a high incidence of brown dwarf 
binaries would be explicable only if the core mass function 
extended to sufficiently low masses (i.e. close to or below the 
hydrogen burning mass limit). Thus the detection of a large 
population of binary brown dwarfs would imply that most 
brown dwarfs arose as the majority component of low mass 
cores, rather than being the low mass ejectae from much 
more massive cores. Clearly, a firmer estimate of the brown 
dwarf binary fraction would place important constraints on 
formation models. 

It is also worth noting that, in these models with min- 
imum core mass of 0.25 M0, the fraction of brown dwarfs 
that remain in bound systems is approximately 25% (75% 
are ejected as singles). This is in clear contrast to objects 
above the hydrogen burning limit, which display a bounded- 
ness fraction close to 90%. 

The kinematics of sub-stellar objects are also likely to 
reflect their origin. Although there is no appreciable depen- 
dence of the final velocity of the singles on mass, brown 
dwarfs acquire as a class higher velocities than stellar ob- 
jects, including binary and multiple systems (sec Figure 14; 
symbols and lines £is in Figure 7). With velocities of a few 
km s~^, it is possible that ejected brown dwarfs would be 
detectable in the vicinity of Class O objects, as suggested by 



Reipurth & Claxke (2001). Relatively few brown dwarfs have 
velocities in excess of 5 km s~^, however, so that the bulk 

of brown dwarfs in the Pleiades would have been retained, 
given that the escape velocity in the young Pleiades is likely 
to have exceeded this value. On the other hand, smaller stel- 
lar associations such as Taurus- Auriga could have already 
lost a significant fraction of its sub-stellar population and 
therefore, would show a relative excess of binary stars. The 
findings of Briceno et al. (2002) axe in line with this argu- 
ment: their survey indicates that Taurus has ~ 2 x fewer 
brown dwarfs at 0.02-0.08 than the Trapezium, whereas 
it is well known (Chez et al. 1993; Leinert et al. 1993; Simon 
et al. 1995, Kohler & Leinert 1998) that the binary frequency 
among T Tauri stars in Taurus is enhanced by a factor of 
two compared to solar type stars on the main sequence. 



8 CONCLUSIONS 

We have modelled the decay of non-hierarchical gas rich sys- 
tems comprising 5 seed protostars in order to study the mass 
function, binarity and kinematics of the resulting popula- 
tion. This study can be seen as complementary to that of 
Bate, Bonnell & Bromm (2002a, b; 2003) who have mod- 
elled similar processes in the context of the calculation of 
the global evolution of a molecular cloud. Although the ap- 
proach used here has the disadvantage that its initial con- 
ditions are artificial, and that it evidently cannot treat any 
possible interactions between star forming cores, the com- 
putational economy effected in not treating the whole cloud 
ensures a) that one can obtain statistically significant popu- 
lations of stars for a given initial condition and hence isolate 
the factors that determine the properties of the resulting 
stars and b) that one can pursue the evolution to comple- 
tion, i.e. to the point where all the gas has been accreted 
and where the system has decayed into a hierarchical mul- 
tiple system. This latter feature allows us, in contrast to 
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the calculations of Bate et al., to explore the properties of 
wider binaries, which remain undecayed in the Bate et al 
simulations. 

Our chief conclusions are as follows: 

(i) The mass function of the resulting stars is primarily 
determined by the mass function of parent cores, being com- 
parable to this core mass function down to stellar masses 
that are similar to the minimum core mass, and declining at 
masses below this. The reason for this is that the division of 
mass for given core mass is essentially bimodal - there is a 
high mass peak corresponding to a binary pair and a lower 
mass peak corresponding to ejected singles or wide tertiary 
companions. When this is convolved with a core mass func- 
tion that declines steeply towards high masses, the majority 
of stars of given mass belong to the high mass peak of low 
mass cores rather than the low mass peak of rare, high mass 
cores. Consequently, the stellar mass function follows the 
core mass function. Thus the observed similarity between 
the core mass function and the stellar IMF (Motte et al. 
1998) cannot of itself be used to disprove the hypothesis that 
most stars arise in small non-hierarchical multiples. 

(ii) Stars much less massive than the minimum core mass 
are unlikely to be binary primaries. This can be understood 
in terms of the evolutionary pattern described above. Thus 
if a large fraction of brown dwarfs turn out to be binary 
primaries, this must imply that the parent core distribu- 
tion extends well into the substcUar regime. In this case, 
although brown dwarfs would form in non-hieraxchical mul- 
tiples, the brown dwarf binaries would not be ejected from 
these systems but would represent the highest mass ob- 
jects in their parent cores. Thus a firm evaluation of the 
binary fraction amongst brown dwarfs will place important 
constraints on the properties of their parent cores. [There 
is the additional possibility that star-disc interactions can 
produce tight brown dwarf-brown dwarf binaries that would 
survive ejection. The lack of discs in the simulations reported 
here do not allow us to model this process. The simulations 
of Bate et al. (which do include discs) indicate that this is 
not common however]. 

(iii) Binary pairs with extreme mass ratios (g < 0.2) are 
not produced in the simulations. More extreme mass ra- 
tios only occur in the case of low mass tertiary companions 
in a wide orbit around a central tight binary. Such a con- 
figuration is the natural consequence of the decay of non- 
hierarchical systems, in which low mass stars may cither be 
ejected or else end up in wide hierarchical orbits. We point 
out that it is an expectation of dynamical decay models 
that where brown dwarfs are found as distant companions 
to normal stars (Gizis et al 2001; Kirkpatrick et al. 2001), 
the central object should normally be a close binary. Over 
all, the distribution of binary mass ratios (Figure 4) is well 
matched to the observations of Mazeh et al (1992) for close 
G dwarf pairs but is probably deficient in low mass compan- 
ions when compared with the statistics of G dwarf binaries 
as a whole (Duquennoy and Mayor 1991). 

(iv) Single stars arc ejected with typical velocities of a 
few km s~^ (Figure 7) irrespective of their mass. The final 
velocities of binary stars with respect to their parent cores 
is typically an order of magnitude lower. In principle, this 
signature of dynamical decay could be sought in star form- 
ing regions, although may be masked by the relative motion 



of the parent cores (typically at a few km s^^). Since the 
binary fraction is lower among low mass stars and brown 
dwarfs than in higher mass stars (see (ii) above), this de- 
pedence of ejection speed on binarity means that in prac- 
tice brown dwarfs are born with a higher velocity dispersion 
(with respect to their parent cores) than higher mass stars. 
With typical velocities of a few km s~^ , the majority of brown 
dwarfs would be expected to be retained in open clusters. 
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